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Abstract 

A  new  method,  called  the  QZ  algorithm,  is  presented  for  the 
solution  of  the  matrix  eigenvalue  problem  Ax  =  XBx  with  general 
square  matrices  A  and  B  .  Particular  attention  is  paid  to  the 
degeneracies  which  result  when  B  is  singular.  No  inversions  of  B 
or  its  submatrices  are  used.  The  algorithm  is  a  generalization  of  the 
QR  algorithm,  and  reduces  to  it  when  B  =  I  .  A  Fortran  program  and 
some  illustrative  examples  are  included. 
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1.  Introduction 


We  shall  be  concerned  with  the  matrix  eigenvalue  problem  of 
determining  the  nontrivial  solutions  of  the  equation 
Ax  =  XBx  , 

where  A  and  B  are  real  matrices  of  order  n  .  When  B  is  nonsingular 


this  problem  is  formally  equivalent  to  the  usual  eigenvalue  problem 


B  ^Ax  =  \x  . 


When  B  is  singular,  however,  such  a  reduction  is  not  possible, 
and  in  fact  the  characteristic  polynomial  det(A-KB)  is  of  degree  less 
than  n  ,  so  that  there  is  not  a  complete  s  of  eigenvalues  for  the 
problem.  In  some  cases  the  missing  eigenvalues  may  be  regarded  as 
"infinite".  In  other  cases  the  entire  problem  may  be  poorly  posed.  The 
term  infinite  eigenvalue  is  justified  by  the  fact  that  if  B  is  perturbed 
slightly  so  that  it  is  no  longer  singular,  there  may  appear  a  number  of 
large  eigenvalues  that  grow  unboundedly  as  the  perturbation  is  reduced  to 
zero.  However,  if  det(A-\B)  vanishes  identically,  say  when  A  and  B 
have  a  common  null  space,  then  any  \  may  be  regarded  as  an  eigenvalue. 
Such  problems  have  unusually  pathological  features,  and  we  refer  to  them 
as  "ill-disposed"  problems. 
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In  numerical  work  the  sharp  distinction  between  singular  and  non¬ 
singular  matrices  is  blurred,  and  the  pathological  features  associated 
with  singular  B  carry  over  to  the  case  of  nearly  singular  B  .  The 
object  of  this  paper  is  to  describe  an  algorithm  for  computing  the 
eigenvalues  and  corresponding  eigenvectors  that  is  unaffected  by  nearly 
singular  B  .  The  algorithm,  the  heart  of  which  we  call  the  QZ-algorit'nm, 
is  essentially  an  iterative  method  for  computing  the  decomposition 
contained  in  the  following  theorem  [10]. 


Theorem.  There  are  unitary  matrices  Q,  and  Z  so  that  QAZ  and  QBZ 

are  both  upper  triangular. 

We  say  that  the  eigenvalue  problems  QAZy  =  AQBZy  and  Ax  =  \Bx  are 
unitarily  equivalent.  The  two  problems  obviously  have  the  same  eigenvalues, 
and  their  eigenvectors  are  related  by  the  equation  x  =  Zy  . 

The  algorithm  proceeds  in  four  stages.  In  the  first,  which  is  a 
generalization  of  the  Householder  reduction  of  a  single  matrix  to 
Hessenberg  form  [4,5],  A  is  reduced  to  upper  Hessenberg  form  and  at  the 
same  time  B  i<.  reduced  to  upper  triangular  form.  In  the  second  step, 
which  is  a  generalization  of  the  Francis  implicit  double  shift  QJR  algorithm 
[3>8],  A  is  reduced  to  quasi -triangular  form  while  the  triangular  form 
of  B  is  maintained.  In  the  third  stage  the  quasi-triangular  matrix  is 
effectively  reduced  to  triangular  form  and  the  eigenvalues  extracted.  I:. 
the  fourth  stage  rhe  eigenvectors  are  obJ  „ined  from  the  triangular  matrices 
and  then  transformed  back  into  the  original  coordinate  system. 

The  transformations  used  in  reducing  A  and  B  are  applied  in  such 
a  way  that  Wilkinson's  general  analysis  of  the  roundoff  errors  in  unitary 
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transformations  f 11 ]  shows  that  the  computed  matrices  are  exactly  unitarily 
equivalent  to  slightly  perturbed  matrices  A+E  and  B+F  .  This  means 
that  the  computed  eigenvalues,  which  are  the  ratios  of  the  diagonal  elements 
of  the  final  matrices,  are  the  exact  eigenvalues  of  the  perturbed  problem 
(A+E)x  =  \(B+F)x  .  If  an  eigenvalue  is  well  conditioned  in  the  sense  that 
it  is  insensitive  to  small  perturbations  in  A  and  B  (see  [10]  for  a 
detailed  analysis),  then  it  will  be  computed  accurately.  This  accuracy 
is  independent  of  the  singularity  or  nonsingularity  of  B  . 

The  use  of  unitary  transformations  in  the  reduction  also  simplifies 
the  problem  of  convergence:  a  quantity  may  be  set  to  zero  if  a  perturrat icn 
of  the  same  size  can  be  tolerated  in  the  original  matrix. 

Our  computer  program  does  not  actually  produce  the  eigenvalues  <V 
but  instead  returns  Oh  and  ,  the  diagonal  elements  of  the  triangular 
matrices  QAZ  and  QBZ  .  The  divisions  in  L  -  ay/p^  beCome  "the 
responsibility  of  the  program’s  user.  We  emphasize  this  point  because 
the  ay  and  contain  more  information  than  the  eigenvalues  themselves. 

Since  our  algorithm  is  an  extension  of  the  QR  algorithm,  the  well 
known  properties  of  the  QR  algorithm  apply  to  describe  the  behavior  of 
our  algorithm. 

In  their  survey  article  [9],  Peters  and  Wilkinson  describe  another 
approach  for  the  case  when  B  is  nearly  singular.  In  their  method  on^ 
computes  an  approximate  null  space  for  B  and  removes  it  from  the  problem. 
The  tecnnique  is  reapplied  to  the  deflated  problem,  and  so  on  until  a  well 
conditioned  problem  is  obtained.  The  method  has  the  crucial  drawback  that 
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one  must  determine  the  rank  of  B  .  If  a  wrong  decision  is  reached, 
the  well-conditioned  eigenvalues  may  be  seriously  affected. 

The  special  case  where  A  is  symmetric  and  B  is  positive  definite 
has  been  extensively  treated.  For  the  case  of  well-conditioned  B  the 
"Cholesky -Wilkinson"  method  [ 6  ]  enjoys  a  well  deserved  popularity. 

A  modification  of  this  algorithm  for  band  matrices  is  given  by  Crawford  [  1 
A  variant  of  the  Peters-Wilkinson  method  for  nearly  semidefinite  B  has 
been  given  by  Fix  and  Heiberger  [2].  Although  our  method  does  not 
preserve  symmetry  and  is  consequently  more  time  consuming  than  these 
algorithms,  its  stability  may  make  it  preferable  when  B  is  nearly 
semidefinite. 
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2.  Reduction  to  Hessenberg-Triangular  Form 

In  this  section  we  shall  give  an  algorithm  whereby  A  is  reduced 
to  upper  Hessenberg  form  and  simultaneously  B  is  reduced  to  triangular 
form.  While  a  treatment  of  the  reductions  in  this  and  the  following  sections 
can  be  given  in  terms  of  standard  plane  rotations  and  elementary 
Hermitian  matrices,  we  find  it  convenient  from  a  computational  point  of 
view  to  work  exclusively  with  a  modified  form  of  the  elementary  Hermit ians. 
Accordingly,  we  introduce  the  following  notation. 

By  V ( k)  we  mean  the  class  of  symmetric,  orthogonal  matrices  of 
the  form 


_  .  J- 

I  +  vu 
T 

where  v  u  =  -2  ,  v  is  a  scalar  multiple  of  u  ,  only  components 

k, k+1, . . .,k+r-l  of  u  are  nonzero,  and  vy  =  1  .  Given  any  vector  x  , 

it  is  easy  to  choose  a  member  Q  of  Vr(k)  so  that 

T 

Qx  -  x  +  (u  x)v 

has  its  k+i, . .  .,k+r-l  components  equal  to  zero,  its  k-th  component  changed 
and  all  other  components  unchanged.  Since  u^  =  1  ,  the  computation  of 
Qy  for  any  y  requires  only  2r-i  multiplications  and  2r-l  additions. 
(In  particular,  use  of  a  matrix  in  V2  requires  only  5  multiplications 
instead  of  the  U  required  by  a  standard  plane  rotation.) 

••’or  the  moat  part,  we  shall  use  only  matrices  in  V2  and  y_  . 

When  a  matrix  Q  in  9/^(k)  premultipliss  a  matrix  A  only  rows 
k  ,  k+1  ,  and  k+2  in  QA  are  changed.  If  the  elements  k  ,  k+1  , 
and  k+2  in  a  col.unn  of  A  are  zero,  they  remain  zero  in  QA  .  Likewise, 
if  2cVj(k )  ,  only  columns  k  ,  k+1  ,  and  k+2  are  changed  in  AZ  .  If 
some  row  has  elements  k  ,  k*-l  ,  and  k*2  zero,  then  uhey  remain  zero 
in  AZ  .  Similar  considerations  hold  for  the  class  y_  . 


k 
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All  our  transformations  will  be  denoted  by  Q*s  and  Z' s  with 
various  subscripts.  The  Q’s  will  always  be  premultpliers,  that  is 
row  operations.  The  Z's  will  always  be  postmultipliers,  or  column 
operations.  The  letter  Q  is  being  used  in  its  traditional  role  to 
denote  orthogonal  matrices.  The  letter  Z  was  chosen  to  denote  orthogonal 
matrices  which  introduce  zeros  in  strategic  locations. 

The  first  step  in  the  reduction  is  to  reduce  B  to  upper  triangular 
form  by  T'r emult iplic at ion  by  Householder  reflections.  The  details  of 
this  reduction  are  well  known  (e.g.  see  [^,11])  and  we  confine  ourselves 
to  a  brief  description  to  illustrate  our  notation.  At  the  k-th  stage  of 
the  reduction  (illustrated  below  for  k  -  5  and  n  =  5  ) >  the  elements 
below  the  first  k-1  diagonal  elements  of  B  are  zero. 

X  X  X  X  X 

0  x  x  x  x 

0  0  x  x  x 

0  0  x1  x  x 

0  0  x1  x  x 

Each  x  represents  an  arbitrary  nonzero  element.  Each  x^  represents 
an  element  to  be  annihilated  in  the  next  step.  A  matrix 
is  chosen  to  annihilate  b  and  B  is  overwritten 

hy  Q,  B  ,  giving  a  matrix  of  the  form  illustrated  below. 

iC 

X  X  X  X  X 

0  x  x  x  x 

0  0  x  x  x 

0  0  0  x  x 

0  0  0  x1  x 

This  process  is  repeated  until  k  =  n-l  .  Of  course  A  is  overwritten 

by  VlVs'  '*Q1A  * 
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After  this  reduction,  A  and  B  have  the  form 


A 


B 


X 

X 

X 

X 

X 

X 

X 

X 

X 

V 

X 

X 

X 

X 

X 

0 

X 

X 

X 

X 

X 

X 

X 

X 

X 

0 

0 

X 

X 

V 

X 

X 

X 

X 

X 

0 

0 

0 

X 

X 

1 

X 

X 

X 

X 

X 

0 

0 

0 

0 

X 

The  problem  now  is  reduce  A  to  upper  Hessenberg  form  while 
preserving  the  triangularity  of  B  .  This  is  done  as  follows  ( for  k  =  5  ) • 


First  QeV2(^) 

is  determined  to 

annihilate  a 

51  • 

The 

matrices 

QA 

and 

QB  ,,  which  overwrite 

A 

and  B 

then  have  the 

form 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

0 

X 

X 

X 

X 

X 

X 

X 

X 

X 

0 

0 

X 

X 

X 

X 

X 

X 

X 

X 

0 

0 

0 

X 

X 

0 

X 

X 

X 

X 

0 

0 

0 

1 

x" 

V 

The  transformation  has  introduced  a  nonzero  element  on  the  ( 5, -position 
of  B  .  However,  a  ZeW^k)  can  be  used  to  restore  the  zero  without 
disturbing  the  zero  introduced  in  A  . 

This  step  is  typical  of  all  the  others.  The  elements  of  A  are 
annihilated  by  Q's  in  the  order  illustrated  below. 


X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

x5 

X 

X 

X 

X 

2 

5 

X 

X' 

X 

X 

X 

1 

h 

6 

X 

X 

X 

X 

X 
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As  each  element  of  A  is  annihilated,  it  introduces  a  nonaero  element 
on  the  subdiagonal  of  B  ,  »hich  is  immediately  annihilated  by  a  suitably 

chosen  z  .  The  entire  algorithm,  including  the  Householder  triangulari- 
zation  of  B  may  be  summed  up  as  follows. 

1)  For  k  =  1,2, . . .,n-l 

1)  Choose  5  fj( n-kfl(k)  to  annihilate  b,  _  ,  ,b 

K  n  K+i  k+1,  k  k+2,  k'  ' -; 

2)  B  -  QkB  ,  A  -  QkA 

2)  For  k  =  1,2,  ...,n-2 

1)  For  t  =  n-l,n-2,  ...,k+l 

1)  Choose  Q  eY  (i)  to  annihilate  a.  , 

d  i+l,k 

2)  A  -  QfcfA  ,  B  -  QWB 

3)  Choose  Z ktW2(£)  to  annihilate  b£+1  £ 

The  C°mplete  reduCtion  re<Juires  about  ^  n3  multiplications, 

^  "3  additions  Md  «2  square  roots.  If  eigenvectors  are  also  to  be 
computed,  the  product  of  the  Z’s  must  be  accumulated.  This  requires 
an  additional  §  n3  multiplications  and  §  n3  additions.  The  product 
of  the  Q's  is  not  required  for  the  computation  of  eigenvectors. 
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In  this  and  the  next  section  we  assume  that  A  is  upper  Hessenberg 
and  B  is  upper  triangular.  In  this  section  we  shall  propose  an  iterative 
technique  for  reducing  A  to  upper  triangular  form  while  maintaining  the 
triangularity  of  B  .  The  idea  of  our  approach  is  to  pretend  that  B  is 
nonsingular  and  examine  the  standard  QR  algorithm  for  C  =  AB  1  .  The 
manipulations  are  then  interpreted  as  unitary  equivalences  on  A  and  3  . 

Specifically  suppose  that  one  step  of  the  QR  algorithm  with  'shift  k 
is  applied  to  C  .  Then  Q  is  determined  as  an  orthogonal  transformation 
such  that  the  matrix 
(3.1)  R  =  Q(C-kl) 

is  upper  triangular.  The  next  iterate  C*  is  defined  as 
C‘  =  RQT  +  kl  =  QCQT  . 

If  we  set 

A’  =  QAZ 


and 


B’  =  QBZ  , 

where  Z  is  any  unitary  matrix,  then 

A'B'"1  =  QAZZTB"V  =  QAB'V  =  C' 


The  matrix  Q  is  determined  by  the  requirement  that  R  be  upper 
triangular.  We  choose  Z  so  that  A’  is  upper  Hessenberg  and  3'  is 
upper  triangular.  This  insures  that  the  nice  distribution  of  zeros, 
introduced  by  the  algorithm  of  Section  2,  is  preserved  by  the  QZ  step. 
Thus  a  tentative  form  of  our  algorithm  might  read 


Q 
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1)  Determine  Q  so  that  QC  is  upper  triangular, 

2)  Determine  Z  so  that  QAZ  is  upper  Hessenberg  and  QBZ  is  upper 
triangular , 

3)  A  -  QAZ  ,  B  -  QBZ  . 

The  problem  is  then  to  give  algorithms  for  computing  Q  and  Z  which  do 
not  explicitly  require  C  =  AB_1  . 

The  determination  Q  is  relatively  easy.  For  from  (3.1)  and  the 
definition  of  C  it  follows  that 


(3.2) 


Q(A-kB)  =  RB  =  S  . 


Since  R  and  B  are  upper  triangular,  so  is  S  .  Thus  Q  is  the  unitary 
matrix  that  reduces  A-kB  to  upper  triangular  form.  Since  A-kB  is 
upper  Hessenberg,  Q  can  be  expressed  in  the  form 


(3.3) 


^•WW-^I  ' 


where  Q^V^fk)  . 

To  calculate  Z  we  apply  Q  in  its  factored  form  (3. 3)  to  B  and 
determine  Z  in  a  factored  form  so  that  B  stays  upper  triangular. 


QjB 

has 

the  form 

(k  = 

5) 

X 

X 

X 

X 

X 

1 

X 

X 

X 

X 

X 

0 

0 

X 

X 

X 

0 

0 

0 

X 

X 

0 

0 

0 

0 

X 

If  QjB  is  postmult iplied  by  a  suitable  Z^e?/2(l)  the  nonzero  element 
below  the  diagonal  can  be  removed.  Similarly  QgQ-^BZ.^  has  the  form 


a 

i 


X 


0 

0 

0 

0 


x 

X 

x' 

0 

0 


X 

X 

X 

0 

0 


X  X 
X  X 
X  X 
X  X 

0  x 


and  the  offending  nonzero  element  can  be  removed  by  a  . 

Proceeding  in  this  way,  we  construct  Z  in  the  form 

z  =  z1zg...zh.1  , 

where  . 

Although  QBZ  is  upper  triangular,  it  is  not  at  all  clear  that 
QAZ  is  upper  Hessenberg.  To  see  that  it  is,  rewrite  equation  (3*2) 
in  the  form 

(3 -JO  QAZ  =  SZ  +  kQBZ  . 

From  the  particular  form  of  Z  and  the  fact  that  S  is  upper  triangular, 
it  follows  that  SZ  is  upper  Hessenberg.  Thus  (5*^)  expresses  QAZ  as 
the  sum  of  an  upper  Hessenberg  and  an  upper  triangular  matrix.  In 
fact  (3.1*)  represents  a  computationally  convenient  form  for  computing  QAZ  • 
We  summarize  as  follows. 

1)  Determine  Q  =  *-^i  (\eV2(k))  so  that 

S  =  Q(A-kB)  is  upper  triangular. 

2)  Determine  Z  =  ^i^2*”^n-l  so  that  3f  =  QBZ 

is  upper  triangular 

5}  A:  =  SZ  +  kB* 


If  this  algorithm  is  applied  iteratively  with  shifts  k^,k^, ...  , 
there  result  sequences  of  matrices  A^,A^y  ...  ,  B-^Bg, •••  ,  and 
C1,C2' ***  satisfying 
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Vl  -  ^vZ 


Vl  =  ^ 


B  ,  =QB  Z 
\H-1  Hv 


Cv  =  AvV 


provided  B^  is  nonsingular.  The  matrices  Ay  are  upper  Hessenberg 

and  the  B^  are  upper  triangular.  The  Cv  are  the  upper  Hessenberg 

matrices  that  would  result  from  applying  the  QR  algorithm  with  shifts 

k-^kg, .  .  to  .  As  Cv  tends  to  upper  triangular  form,  so  must  Ay  , 

since  B  is  upper  triangular. 

Most  of  the  properties  of  the  QR  algorithm  carry  over  to  the  QZ 

algorithm.  The  eigenvalues  will  tend  to  appear  in  descending  order  as 

(v} 

one  proceeds  along  the  diagonal.  The  convergence  of  a^  ^  1  to  zero 

may  be  accelerated  by  employing  one  of  the  conventional  shifting  strategies. 

(v) 

Once  a^  ^  ^  becomes  negligible  one  can  deflate  the  problem  by  working 

with  the  leading  principal  submatrices  of  order  n-1  .  If  sane  other 

( v} 

subdiagonal  element  of  A  ,  say  a,  ,  becomes  negligible,  one  can 
effect  a  further  savings  by  working  with  rows  and  columns  l  through  n  . 
Because  we  have  used  unita'y  transformations,  an  element  of  A  or  B. 

V  v> 

can  be  regarded  as  negligible  if  a  perturbation  of  the  same  size  as  the 
element  can  be  tolerated  in  A^  or  B1  . 

The  algorithm  given  above  is  potentially  unstable.  If  k  is  large 
compared  with  A  and  B  ,  the  formula  (3.k)  will  involve  subtractive 
cancellation  and  A'  will  be  computed  inaccurately.  Since  the  shift  .< 
approximates  the  eigenvalue  currently  being  found  and  the  problem  may 
have  very  large  eigenvalues,  there  is  a  real  possibility  of  encountering 
a  large  shift.  Fortunately  the  large  eigenvalues  tend  to  be  found  last 
so  that  by  the  time  a  large  shift  emerges  the  small  eigenvalues  will  have 
been  computed  stably.  (The  large  eigenvalues  are  of  course  ill-conditioned 
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and  cannot  be  computed  accurately.)  To  be  safe  one  might  perform  the 
first  few  iterations  with  a  ?  ro  shift  in  order  to  give  the  larger 
eigenvalues  a  chance  to  percolate  to  the  top. 
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4.  Implicit  Shifts 


The  potential  instability  in  the  explicit  algorithm  results  from 
the  fact  that  we  have  used  formula  (3*4)  rather  than  unitary  equivalences 
to  compute  A’  .  One  way  out  of  this  difficulty  is  to  generalize  the 
implicit  shift  method  for  the  Q,R  algorithm  to  the  QZ  algorithm  so  that 
both  A'  find  B'  are  computed  by  unitary  equivalences.  The  implicit 
shift  technique  has  the  additional  advantage  that  it  can  be  adapted  to 
perform  two  shifts  at  a  time.  For  real  matrices  this  means  that  a  double 
shift  in  which  the  shifts  are  conjugate  pairs  can  be  performed  in  real 
arithmetic. 

Since  we  are  primal’ ily  interested  in  real  matrices,  we  will  concentrate 
on  double  shifts.  The  method  is  based  on  the  following  observation. 

Suppose  that  A  is  upper  Hessenberg  and  B  is  upper  triangular  and 
nonsingular.  Then  if  Q  and  Z  are  unitary  matrices  such  that  Q,AZ 
is  upper  Hessenberg  find  QBZ  is  upper  triangular,  then  Q,  is  determined 
by  its  first  row.  In  fact,  AB  and  QAB  Q  are  both  upper  Hessenberg, 
so  that,  by  the  theorem  on  page  352  of  [11],  Q,  is  determined  by  its 
first  row. 

Thus  we  must  do  two  things.  First,  find  the  first  row  of  3.  . 

Second,  determine  Q,  and  Z  so  that  Q,  has  the  correct  first  row, 

QAZ  is  upper  Hessenberg,  and  QBZ  is  upper  triangular.  The  first  part 
is  relatively  easy.  The  first  row  of  Q  is  the  first  row  that  would  be 
obtained  from  a  double  shifted  QR  applied  to  AB-'*'  .  Since  A  is 
upper  Hessenberg  and  B  upper  triangular,  it  is  easy  to  calculate  the 
first  two  columns  of  AB  ~  .  But  these,  along  with  the  shifts,  completely 
determine  the  first  row  of  Q,  .  Only  nonsingularity  of  the  upper  2-by-2 
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submatrix  of  B  is  actually  required  here.  If  either  b1  -|  or  b^ 
is  too  small,  so  that  this  submatrix  is  nearly  singular,  a  type  of 
deflation  can  be  carried  out.  We  will  return  to  this  point  later. 

The  second  part  is  a  little  more  difficult,  and  is  really  the  crux 
of  the  algorithm  since  it  retains  the  Hessenberg  and  triangular  forms. 
Only  the  first  three  elements  of  the  first  row  of  Q,  are  nonzero.  Thus, 


if  Qx 

is 

a  matrix 

in 

tf3(l) 

with  the 

same 

first 

row 

of 

and  (^B 

have  the  following  form  (when 

n  = 

6  ) 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

2 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

1 

X 

1 

X 

X 

X 

X 

X 

0 

0 

X 

X 

X 

X 

C 

0 

0 

X 

X 

X 

0 

0 

0 

X 

X 

X 

0 

0 

0 

0 

X 

X 

0 

0 

0 

0 

X 

X 

0 

0 

0 

0 

0 

X 

As  in  the  standard  implicit  shift  QR  algorithm,  it  is  convenient  to  think 
of  as  the  reflection  which  annihilates  two  of  the  three  nonzero 
elements  in  a  fictitious  "zeroth"  column  of  A  . 

We  must  reduce  Q^A  to  upper  Hessenberg  and  Q^B  to  upper  triangular 
by  unitary  equivalences.  However,  we  may  not  premultiply  by  anything  which 
affects  the  first  row.  This  is  done  as  follows.  The  matrix  Q^B  has 
three  nonzero  elements  outside  the  triangle.  These  can  be  annihilated 
by  two  Z*s  ,  a  Zj,  in  3^(1)  which  annihilates  the  (5,1)  and  (5,2) 
elements  and  then  a  Z£  which  annihilates  the  resulting  (2,1)  element. 
Let  Z^  =  Zj_Z!^  .  Then  Q-|BZ^  is  upper  triangular.  Applying  Z^  to 
0  -  A  gives  0  .  AZ  +*  n  +  V>  o  -Pr»1  1  rnH  n  rr  *Pnr*m 

*'2  *  O'*  ^"*1  ',**w 
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X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

1 

X 

X 

X 

X 

X 

X 

1 

X 

X 

X 

X 

X 

X 

0 

0 

0 

X 

X 

X 

0 

0 

0 

0 

X 

X 

This  is  multiplied  by 

Q2  in 

V;(2) 

that 

annihilates  the 

(5,1) 

and 

(4, l)  elements. 

Then 

and 

WZ1 

have  the  forms 

X 

X 

X 

X  X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X  X 

X 

0 

X 

X 

X 

X 

X 

0 

X 

X 

X  X 

X 

0 

X2 

X 

X 

X 

X 

0 

X 

X 

X  X 

X 

0 

X^ 

1 

X 

X 

X 

X 

0 

0 

0 

X  X 

X 

0 

0 

0 

0 

> 

X 

0 

0 

0 

0  x 

X 

0 

0 

0 

0 

0 

X 

The  first  columns  are  now  in  the  desired  form.  The  nonzero  elements 
outside  the  desired  structure  have  been  "chased"  into  the  lower  5-by~5 
submatrices . 

Now,  postmult iply  by  ,  a  product  of  a  matrix  in  ^(2)  and  a 
matrix  in  $^(2)  that  reduces  the  current  B  to  triangular  form.  Then 
premultiply  by  in  5^(3)  to  annihilate  two  elements  outside  the 

Hessenberg  structure  of  the  resulting  A  . 

The  process  continues  in  a  similar  way,  chasing  the  unwanted  nonzero 
elements  towards  the  lower,  right-hand  comers.  It  ends  with  a  slightly 
simpler  step  which  uses  Qn  ^  in  V2(n-l)  to  annihilate  the  (n,n-2) 
element  of  the  current  A  ,  thereby  producing  a  Hessenberg  matrix,  and 
Zn_2  in  ^(n-l)  which  annihilates  the  (n,n-l)  element  of  the 
current  B  ,  producing  a  triangular  B  but  not  destroying  the  Hessenberg  A  . 


1 6 
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We  are  now  in  a  position  to  summarize  the  double  implicit  shift 
method.  It  is  understood  that  A  and  B  are  to  be  overwritten  by 
the  transformed  matrices  as  they  are  generated. 


1)  Compute  a1Q}  a^,  and  a^  by  (t.l) . 

2)  For  k  =  1,2, . ..,n-2 

a)  Determine  Q^?/^(k)  to  annihilate  a^^  ^  ^  and  a^  ,  k  ^ 

b)  Determine  Z^ejf^(k)  to  annihilate  b^+0  k+1  and  b}4^  ^  . 

c)  Determine  Z!^U^(k)  to  annihilate  ^  . 

3)  Determine  Q  ^eV-fn-l)  to  annihilate  a  „  . 

n-1  2V  '  n,n-2 

b)  Determine  Z  .,e?/,-(n-l)  to  annihilate  b  . 

'  n-1  2V  '  n,n-l 


For  each  k  ,  determination  of  Q,^  requires  a  few  multiplications 


and  one  square  root.  Application  of  to  both  A  and  B  requires 

about  10(n-k)  multiplications.  The  work  involved  with  each  Z!  is 


the  same.  Application  of  Z£  requires  only  about  6(r.-k)  multiplications. 


The  number  of  additions  is  about  the  same.  Summing  these  for  k  from 

2  2 

1  to  n-1  gives  a  total  of  about  13n  multiplications,  15n  additions 


and  3n  square  roots  per  double  iteration. 


By  way  of  comparison,  for  the  double  shift  QR  algorithm  as  implemented 

T 

in  "  hqr  ",  Z £  becomes  simply  and  Z £  is  not  used.  Furthermore,  the 

transformations  are  carried  out  on  only  one  matrix.  Consequently,  each 

o  2 

double  iteration  requires  about  5n~  multiplications,  5n  additions  and 


n  square  roots.  Thus  the  QZ  algorithm  applied  on  two  matrices  can  be 
expected  to  require  roughly  2.6  times  as  much  work  per  iteration  as 
the  QR  algorithm  on  a  single  matrix. 


r,AV«V'>.,  A,  t.,'w"SS', 
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In  order  to  obtain  eigenvector  the  Q's  are  ignored  and  the  Z's 

O  O 

accumulated.  This  requires  about  8n  more  multiplications  and  8n 
more  additions  per  double  iteration. 

There  is  one  difficulty.  The  formulas  for  a^Q  ,  ,  and 

are  not  defined  when  b^  and  b^  are  zero.  Moreover,  as  b^  and 
hgg  approach  zero  the  terms  that  determine  the  shift  (terms  involving 
ann  ’  fcnn  ’  ei:c*)  become  negligible  compared  to  the  other  terms,  so  that 
the  effect  of  the  shift  is  felt  only  weakly. 

Part  of  the  solution  to  this  difficulty  is  to  deflate  from  the  top. 
If  b11  is  negligible  it  may  be  set  to  zero  to  give  the  forms  for  A 
and  B  (n  =  4) 


X 

X 

X 

X 

0 

X 

X 

X 

X 

X 

X 

X 

0 

X 

X 

X 

0 

X 

X 

X 

0 

0 

X 

X 

C 

0 

X 

X 

0 

0 

0 

X 

A  Q,  in  V2(l)  can  then  be  used  to  annihilate  the  (2,1)  element  of  A  , 
which  deflates  the  problem. 

The  rest  of  the  solution  lies  in  recognizing  that  thei-e  .is  not  muo: 
of  a  problem.  If  b.^  and  b^  are  small  then  the  problem  has  large 
eigenvalues.  We  have  already  observed  that  the  larger  eigenvalues  f end 
to  emerge  at  the  upper  left,  and  the  larger  the  eigenvalue,  the  swifter 
its  emergence.  Moreover  the  speed  will  not  be  affected  by  a  small  shift. 
This  means  that  whenever  the  implicit  shift  is  diluted  by  a  small  b(, 
or  >  the  algorithm  is  none  the  less  profitably  em]J.ove<?  in  finding 

a  large  eigenvalue. 
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5-  Furthe.  1  eduction  of  the  Quasi-Triangular  Form 

The  result  of  the  algorithm  described  so  far  is  in  an  upper 
triangular  matrix  B  and  a  quasi-upper  triangular  matrix  A  in  which 
no  two  consecutive  subdiagonal  elements  are  nonzero.  This  means  that 
the  original  problem  decomposes  into  one  by  one  and  two  by  two  subproblems. 

The  eigenvalues  of  the  one  by  one  problems  are  the  ratios  of  the  corres¬ 
ponding  diagonal  elements  of  A  and  B  .  The  eigenvalues  of  the  two  by 
two  problems  might  be  calculated  as  the  roots  of  a  quadratic  equation, 
and  may  be  complex  even  for  real  A  and  B  . 

There  are  two  good  reasons  for  not  using  the  quadratic  directly, 
but  instead  reducing  the  two  by  two  problems.  First,  when  A  and  B 
are  real,  the  calculation  of  eigenvectors  is  greatly  facilitated  if  all 
the  real  eigenvalues  are  contained  in  one  by  one  problems.  A  more 
important  second  reason  is  that  the  one  by  one  problems  contain  more 
informat i-»n  then  the  eigenvalues  alone.  For  example,  if  a^  and 
are  small  then  the  eigenvalue  =  a^/b_^  is  ill  conditioned,  however 
reasonable  it  may  appear.  This  reason  obviously  applies  to  complex 
eigenvalues  as  well  as  real  ones.  Accordingly,  we  recommend  that  the 
two  by  two  problems  be  reduced  to  one  by  one  problems  and  that  the 
diagonal  elements,  rather  than  the  eigenvalues,  be  reported. 

Without  loss  of  generality  we  may  consider  the  problem  of  reducin' 
two  by  two  matrices  A  and  B  simultaneously  to  upper  triangular  form 
by  unitary  equivalences.  For  our  purposes  we  may  assume  that  B  is 
upper  triangular. 

Two  special  cases  may  be  disposed  of  immediately.  If  b  is  zero, 

then  a  Qc?/0(l)  may  be  chosen  to  reduce  a to  zero.  The  zero  elements  | 

d  dl 
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of  QB  are  not  disturbed.  Similarly,  if  is  zero,  a  ZGVg(l) 

may  be  chosen  to  reduce  to  zero  without  disturbing  b^  . 

In  the  general  two-by-two  case,  it  is  not  difficult  to  write  down 
formulas  for  the  elements  of  A*  =  QAZ  and  B*  =  QBZ  for  any  Q  and  Z  . 
Moreover,  these  formulas  can  be  arranged  so  that  numerically  one  of 
or  b^  is  effectively  zero.  It  is  not  obvious,  however,  that  the  other 
element  is  numerically  zero,  and  the  effect  of  assuming  that  it  is  be¬ 
setting  it  to  zero  could  be  disastrous.  Consequently,  we  must  consider 
a  somewhat  more  complicated  procedure. 

The  theoretical  procedure  for  reducing  A  to  triangular  form  may  be 
described  as  follows.  Let  X  be  an  eigenvalue  of  the  problem  and  form 
the  matrix  E  =  A-XB  .  Choose  a  Zej^l)  to  annihilate  either  e^ 
or  .  Since  the  rows  of  E  are  parallel,  it  follows  that  whichever 

of  e^  or  egl  is  annihilated  the  other  must  also  be  annihilated. 

Now  choose  QeY2(l)  so  that  either  QAZ  or  QBZ  is  upper  triangular. 

Since  the  first  column  of  QEZ  is  zero  and  QEZ  =  QAZ  -  XQBZ  ,  it  follows 
that,  however  Q  is  chosen,  both  QAZ  and  QBZ  must  be  upper  triangular. 

In  the  presence  of  rounding  error  the  method  of  computing  X  and 
the  choice  of  Z  and  Q  are  critical  to  the  stability  of  the  process. 

A  rigorous  rounding  error  analysis  will  show  that,  under  a  reasonable 
assumption  concerning  the  computed  X  ,  the  process  described  below  is 
stable.  However,  to  avoid  excessive  detail,  we  only  outline  the  analysis. 

We  assume  that  all  computations  are  done  in  floating  point  arithmetic  with 
t  base  p  digits  and  that  the  problem  has  been  so  scaled  that  underflows 
and  over  flows  do  not  occur.  We  further  assume  that  is  not  negligible 

in  the  sense  that  |a21l  <  3_^||A||  >  where  ||*||  denotes,  say,  the  row 
sot.  norm. 


1 

\ 

-! 


/ 
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The  algorithm  for  computing  X  amounts  to  making  an  appropriate 


origin  shift  and  computing  an  eigenvalue  from  the  characteristic  equation. 


It  goes  as  follows 
a 


V-  = 


11 


Jn 


a12  _  al2  ^b12 


22 


=  a, 


22  ^22 


(jl>r 


if  ^22  ° 
=  5Vb22  '  b 


12a21 


llb22 


q  = 


a2iai2 


bnb22 


(5.1) 


=  p  +  q 

=  p.  +  p  +  sign(p)  «/r 


(complex  if  r  <  0  ) 


We  must  now  assume  that  the  computed  k  satisfies  the  equation 
det(A*  -\B»)  =  0  , 

where  ||A-A'  j|  <  ^a||aJ|  and  J|B~B* ((  <  0g|jB||  with  and  smail 
constants  of  order  0  b  .  Define 
E*  =  A*  -\B» 

and  let  E  denote  the  computed  value 
E  ---  fl(A-XB)  . 

Then 


E»  =  E  +  H 


with  jjK|i  <  o  ma>- {jjAjj ,  jkj  jjBjj }  with  o  of  order  (3  b  . 


We  claim  that,  approximately, 


(‘i.2)  jlEji  >  p”b  max{llA||,  M  11B1I) 


22 


I 
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First  we  note  that 


(5-3)  ||E||  >  je  |  =  ja  |  >  p‘ 


by  the  assumption  that  is  significant.  New  assume  that 

||E|j  <  0  ^|Xj  ||B||  .  Then  subtractive  cancellation  must  occur  in  the 
computation  of  e^  ,  e^  ,  and  e22  .  Thus  «  Xb^  ,  a^2  ~  XbJ2 

and  a22  Xb22  .  Hence  we  have  ||Aj|  >  |X|  ||b|{  ,  and,  from  (5*3), 

jjEjj  >  0  ^  |X|  |jBjj  ,  a  contraction. 


0  =  det(E')  =  det(E)  +  (e^h^h^  -  ( e^+h^) h21  +  h^e^  -  h12e21  , . 

Hence 

|det(E)  |  <  pJjEjl  max{||A|J,  |X|  ||B||  }  +  p2[maxf||A||,  |X|  ||b||  }f 

where  p^  and  p2  are. of  order  0  ^  .  From  (5*2)  it  then  follows  that 

|det(E)  |  <  pl|E||  max{||A||,  jX|  ||b|'  } 

“"t 

where  p  is  of  order  0 

Now  consider  the  determination  of  Z  .  Assume  that  the  second  row 
of  E  is  larger  than  the  first.  Then  ZtV2(l)  is  chosen  to  annihilate  e21 
Let  F  =  EZ  Then  f^  is  essentially  zero.  Furthermore,  since  Z  is 
unitary 

1  1  =  |&et(E)  |  <  p||  E |[  max{||A||,  |x| ||b||} 

But  |f22[  ^  ll^ll  since  e2  was  assumed  to  be  the  larger  row, 

||e2l!  =  ||E||  .  Hence  we  have  approximately 

I  .p  |  ^  „ _ fit.  II  )il  11—in 

l-ii  I  ^  P  ■■k»aI||M||,  HlPIIJ  • 


To  choose  Q  ,  let 


C  =  AZ  ,  D  =  BZ  , 
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and  let  f ^  ,  and  d^  be  the  first  columns  of  F  ,  C  ,  and  D  . 

m 

Let  denote  the  second  row  of  Q,  .  If  j|A|j  >  |\  |  ||b||  ,  we  choose  Q 
to  annihilate  d^  .  Numerically  this  means  that 

l4  di!  <  * 

-t  T 

where  a  is  a  constant  on  the  order  of  p  .We  must  show  that  q^c^ 

is  negligible.  But 

\4\\  =  l«j>  fl  +  X  qo  dal 

<  *  IM  Ihe  aj 

<  pmax(!lA!|>|X|||B||)  *  c!MM 

<  (p  +  a)l|A(l  • 

If,  on  the  other  hand,  |M||B|1  >  llA||  ,  we  choose  Q  so  that 

1^2  cl!  <  o-llAjS  . 

It  then  follows  that 


J 

< 

3 

I 

j 

1 


j 


lq2  dll  =  lql  f  "  4  Cll  /  M 

<  pllj"1  max{||A||,|\|||B||}  +  alM“V || 

<  (p  +  a)  IjB'i'i  . 

In  summary,  \  is  computed  using  (5*1),  Z  is  chosen  to  annihilate 
the  first  element  of  the  larger  of  the  two  rows  of  A-M3  and  Q.  is  chosen 
to  annihilate  the  (2,1)  element  of  the  smaller  of  the  two  matrices  A Z 
and  c.BZ  .  In  this  way,  we  can  be  sure  that  the  computed  (2,1)  elements 
of  both  QAZ  and  QBZ  are  negligible. 


In  practice  with  matrices  of  any  ordei,  if  the  transformations  are 
real,  they  are  applied  to  the  entire  matrices.  If  the  transformations  are 


complex,  they  are  used  to  compute  the  diagonal  elements  that  would  result, 
but  are  not  actually  applied.  We  thus  obtain  a  quasi -triangular  problem 
in  which  each  two-by-two  block  is  known  to  correspond  to  a  pair  of  complex 
eigenvalues . 

The  generalized  eigenvectors  of  this  reduced  problem  can  be  found  by 
a  back-substit  .'t-ion  process  which  is  a  straightforward  extension  of  the 
method  used  in  "  hqr2  "  [ 8  ] .  The  vectors  of  the  original  problem  are 
then  found  by  applying  the  tic  emulated  Z's  . 


6.  Some  Numerical  Results 

The  entire  process  described  above  has  been  implemented  in  a  Fortran 
program  [ 7  ]•  There  are  Tour  main  subroutines:  the  initial  reduction  to 
Ressenberg-triangular  form,  the  iteration  itself,  the  computation  of  the 
final  diagonal  elements,  and  the  computation  of  the  eigenvectors.  The 
complete  program  contains  about  600  Fortran  statements,  although  this 
could  be  reduced  somewhat  at  the  expense  of  some  clarity. 

The  numerical  properties  observed  experimentally  are  consistent  with 
the  use  of  unitary  transformations.  The  eigenvalues  are  always  found  to 
whatever  accuracy  is  justified  by  their  condition.  If  an  eigenvalue  and 
eigenvector  are  not  too  "ill-disposed",  then  they  produce  a  small  relative 
residual. 

Similar  numerical  properties  can  not  generally  be  expected  from  any 
algorithm  which  inverts  B  or  any  submatrix  of  B  .  This  is  even  true  of 
2-by-2  submatrices,  as  illustrated  by  the  following  example  due  to 
Wilkinson. 


Here  p  is  about  the  square  root  of  the  machine  precision,  that  is,  u  is 

2 

not  negligible  compared  to  1  ,  but  p  is.  There  is  one  eigenvalue 
near  -2  .  Small  relative  changes  in  the  elements  of  the  matrices  cause 
only  small  relative  changes  in  this  eigenvalue.  The  other  eigenvalue 
becomes  infinite  as  p  approaches  aero.  Great  care  must  be  taken  in 
solving  this  problem  so  that  the  mild  instability  of  the  one  eigenvalue 
does  not  cause  an  inaccurate  result  for  the  other,  stable  eigenvalue. 
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Of  course,  the  use  of  unitary  transformations  makes  our  technique 
somewhat  slower  than  others  which  might  be  considered.  But  the  added 
cost  is  not  very  great.  In  cesting  our  program,  we  solve  problems  of 
order  50  regularly.  A  few  problems  of  orders  greater  than  100  have 
been  run,  but  these  become  somewhat  expensive  when  they  are  merely  tests. 

One  typical  ex&ople  of  order  50  requires  45  seconds  on  Stanford’s 
IBM  560  model  67.  Of  this,  15  seconds  are  spent  in  the  initial  reduction 
29  seconds  are  used  for  the  6l  double  iterations  required,  and  5  seconds 
are  needed  for  the  diagonal  elements  and  eigenvectors.  If  the  eigenvectors 
are  not  needed  and  so  the  transformations  not  saved,  the  total  time  is 
reduced  to  27  seconds.  By  way  of  comparison,  formation  of  B~^A 
a  la  Peters  and  Wilkinson  [9]  and  use  of  Fortran  versions  [12]  of  "orthos" 
[5]  and  "  hqr2  ”  [8]  requires  a  total  of  27  seconds  for  this  example. 

(All  of  these  times  are  for  code  generated  by  the  IBM  Fortran  IV  compiler, 

H  level,  with  the  optimization  parameter  set  to  2  .) 

In  the  examples  we  have  seen  so  far,  the  total  number  of  double 
iterations  required  is  usually  about  1,2  or  1.5  times  the  order  of 
the  matrices.  This  figure  is  fairly  constant,  although  it  is  not  difficult 
to  find  examples  which  require  many  fewer  or  many  more  iterations.  As  a 

rule  of  thumb,  for  a  matrix  of  order  n  the  time  required  on  the  model  -7 

5  5 

is  about  .56  n  milliseconds  if  vectors  are  computed,  .22  n"  milli¬ 
seconds  if  they  are  not. 
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The  example  in  Table  1  is  not  typical,  but  it  does  illustrate 
several  interesting  points.  It  was  generated  by  applying  non -orthogonal 
rank  one  modifications  of  the  identity  to  direct  sums  of  companion  matrices. 
The  companion  matrices  were  chosen  so  that  the  resulting  problem  has 
three  double  roots. 


\  ■  N>  =  " 


x  -  x  _  1  .  /5  , 

~  2  2  1  ’ 

X  _  X  -  1  ft  s 

\  -  \  -  2  '  T  1  • 


The  double  root  at  co  results  from  the  fact  that  B  has  a  double  zero 

eigenvalue.  All  three  roots  are  associated  with  quadratic  elementary 

divisors;  i.e.,  each  root  has  only  one  corresponding  eigenvector.  The 

computed  diagonals  of  the  triangularized  matrices  are  given  in  the  table. 

Note  that  the  four  finite  eigenvalues  are  obtained  with  a  relative  accuracy 
-8 

of  about  10  .  This  is  about  the  square  room  of  the  machine  precision 

and  is  the  expected  behavior  for  eigenvalues  with  quadratic  elementary 

divisors.  The  singularity  of  B  does  not  cause  any  further  deterioration 

in  their  accuracy.  Furthermore,  the  infinite  eigenvalues  are  obtained  from 

the  reciprocals  of  quantities  which  are  roughly  the  square  root  of  the 

machine  precision  times  the  norm  of  B  .  Consequently  we  are  somewhat 

/ 

justified  if  we  claim  to  have  computed  the  square  root  of  infinity.-' 

*T "  “  '  ~  - - 

This  prompts  us  to  recall  the  limericx  which  introduces  George  Gamow’s 
One,  Two,  Three,  Infinity: 

There  was  a  young  fellow  from  Trinity 
Who  tried  To 

But  the  number  of  digits 
Gave  him  such  fidgits 
That  he  gave  up  Math  for  Divinity. 
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mmmm 


m 


50 

-60 

50 

-27 

6 

6 

16 

5 

5 

5 

-6 

5 

38 

-28 

27 

-17 

5 

5 

5 

16 

5 

5 

-6 

5 

27 

-17 

27 

-17 

5 

5 

5 

B  = 

5 

16 

5 

-6 

5 

27 

-28 

38 

-17 

5 

5 

5 

5 

5 

16 

-6 

5 

27 

-28 

27 

-17 

16 

5 

5 

5 

5 

5 

-6 

16 

27 

-28 

27 

-17 

5 

16 

6 

6 

6 

6 

-5 

6 

a  p 

25.768670843143  .2637605112. 10“ 6 

-12.821841071323  . 1312405807. 10"6 

5.814535434181  +  10.071071345641  i  11.629071028730 

5.800765071150  -  10.047220375909  i  11.601530302268 

5.736511506410  +  9.935928843473  i  11.473022854605 

5.510879468089  -  9-5^5322710676  i  11.021758784136 

a/p 

Q 

0.976972281.10 

Q 

-0.976972290.10 

0.49999999310489  +  0.86602543924271  i 
0.49999999310489  -  0.86602543924271  i 
0.50000000689511  +  0.86602536832617  i 
O.5OOOOOOO6895H  -  0.86602536832617  i 


Table  1 
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SUBROUTINE  QZ (ND*N*A*8*EPS» ALFR*ALFI *BETA* ITER* WANTX «X) 

DIMENSION  A (ND  *ND ) *B (ND*NO) ♦ ALFR (N) *ALFI (N) *6ETA(N)*X(ND*ND) 
DIMENSION  ITER <N) 

LOGICAL  WANT X 

A  ANO  B  ARE  N-BY-N  REAL  MATRICES*  STORED  IN  ARRAYS  WITH  NO  ROWS. 
EPS  IS  THE  RELATIVE  PRECISION  OF  ELEMENTS  OF  A  ANO  B. 

FINDS  N  PAIRS  OF  SCALARS*  (ALFA  (M >  *BETMM)  )  SO  THAT 
BETA(M)*a  -  ALFA(M)*B  is  singular. 

The  EIGENVALUES  OF  A*X  -  LAMBDA*B*X  CAN  BE  OBTAINED  BY 
DIVIDING  ALFA(M)  BY  BETA (M) *  EXCEPT  BET  A (M)  MIGHT  BE  ZERO. 

IF  (wANTX)  ALSO  FINDS  CORRESPONDING  EIGENVECTORS. 

USES  ONLY  UNITARY  TRANSFORMATIONS*  NO  INVERSES. 

SO  EITHER  A  OR  B  (OR  BOTH)  MaY  3E  SINGULAR. 

BETA(M)  IS  REAL. 

alfa(m)  is  complex*  real  and  imaginary  parts  in  alfr<m>  anc  alfkm) 

COMPLEX  PAIRS  OCCUR  WITH  ALFA (M ) /BETA (M)  AND  ALFA (M* 1 ) /BETA (M* 1 > 
COMPLEX  CONJUGATES  EVEN  THOUGH  ALFA (M)  AND  ALFA (M* 1 )  ARE  NOT 
NECESSARILY  CONJUGATE. 

USES  ONLY  REAL  ARITHMETIC. 

IF  A  AND  B  WERE  REDUCED  TO  TRIANGULAR  FORM  BY  UNITARY  EQUIVALENCES* 
ALFA  AND  BETA  WOULD  BE  THE  DIAGONALS. 

A  AND  B  ARE  ACTUALLY  REOUCED  ONLY  TO  QUASI “TRIANGULAR  FORM  WITH 
1-BY-T  AND  2-8Y-2  BLOCKS  ON  DIAGONAL  OF  A. 

IF  ALFA(M)  IS  NOT  REAL*  THEN  BETA (M)  IS  NOT  ZERO. 

ITER  IS  TR0U8LE  INDICATOR  AND  ITERATION  COUNTER. 

IF  (ITER(l) .EQ.O)  EVERYTHING  IS  OK. 

ITER (M)  IS  NUMBER  OF  ITERATIONS  NEEDED  FOR  M“TH  EIGENVALUE. 

IF  (ITER(l)  THRU  ITER (M)  .EG.  -l>  THEN  ITERATION  FOR  M-TH 
EIGENVALUE  DID  NOT  CONVERGE  AND  ALFA(1>  THRU  ALFA (M)  AND 
BETA ( 1 )  THRU  8ETA(M)  ARE  PROBABLY  INACCURATE. 

IF  (WANTX)  X ( . *M)  IS  THE  M-TH  REAL  EIGENVECTOR. 

X ( • *M)  AND  X ( • *M* 1 )  ARE  THE  REAL  AND  IMAGINARY  PARTS 
OF  THE  M-TH  COMPLEX  EIGENVECTOR. 

X(.*M)  AND  “X(.*M*D  AND  THE  REAL  AND  IMAGINARY  PARTS 
OF  THE  (M.l)-ST  COMPLEX  EIGENVECTOR. 

VECTORS  NORMALIZED  SO  THAT  LARGEST  COMPONENT  IS  1.  OR  l.*0.I  . 

USES  FOUR  PRIMARY  SUBROUTINES*  QZHES *  QZIT*  QZVAL  AND  QZVEC. 

USES  FOUR  AUX ILL  I ARY  SUBROUTINES*  HSH3*  JSH2*  CHSH2  AND  CDIV. 

USES  TWO  STANDARD  FUNCTIONS*  SORT  AND  ABS. 

AUTHORS*.  C.  B.  MOLER*  STANFORD*  AND  G.  W.  STEWART*  U.  Op  Tf:  XAS 
THIS  VERSION  OATED  7/19/71. 

CALL  OZHES(ND*N*A*B*WANTX*X) 

CALL  QZIT  (ND*N*A*B*EPS*EPSA*EPSB*ITER*WANTX.X) 

CALL  QZVAL  (ND*N*A»B»EPSA»EPSB*ALFR*AI.FI*BETA*WANTX*X) 

IF  (WANTX)  CALL  QZVEC ( NO *N *  A  * d *EPSA .EPSB* ALFR * ALF I ♦ RET A, X ) 

Rt  turn 

(■  hn 
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SUBROUTINE  QZHES <ND*Nf A*B,WANTX*X> 
DIMENSION  A (N0«ND) *B (NOtND) ♦ X (ND fND) 
LOGICAL  WANTX 

INITIALIZE  X*  USEO  TO  SAVE  TRANSFORMATIONS 

IF  (.NOT.WANTX)  GO  TO  10 
DO  3  1=1 »N 
DO  2  J=1 *N 

X(ItJ)  =  0. 

2  CONTINUE 
X(Ifl)  =  1. 

3  CONTINUE 


REDUCE  8  TO  UPPER  TRIANGULAR 

10  NKlsN-1 

DO  100  L=1*NM1 
LI  =  L* 1 

S  =  0. 

DO  20  I=L1 *N 

IF  ( ABS (B ( I »L) ) »GT*S)  S  =  ABS(B(I*L) ) 
20  CONTINUE 

IF  (S.EG.0.1  GO  TO  100 
IF  (ABS(B(L*L)).GT.S)  S  =  ABS (B(L*L) ) 

R  =  0. 

DO  25  I-L*N 

B( I*L>  =  B ( I *L) /S 
R  =  R  ♦  8 < I *L)**2 
25  CONTINUE 

R  =  SQRT(R) 

IF  (P(L.L).LT.O.i  R  =  -R 
BIL.L)  =  B<L»L)  ♦  R 
RrtO  =  R*fi (L  *L ) 

DO  50  J-Ll *N 
T  =  0. 

DO  30  I=L*N 

T  =  T  ♦  B  ( I  *L )  *B  ( I  »  J) 

30  CONTINUE 

T  =  -T/RHO 
DO  40  I=LfN 

B(ItJ)  =  B (I  * J)  ♦  T*B(ItL) 

40  CONTINUE 

50  CONTINUE 

DO  80  J=1.N 
T  =  0. 

DO  60  I=L*N 

T  =  T  ♦  B ( I *L ) *A ( I « J) 

60  CONTINUE 

T  =  -T/RhG 
DO  70  I  =L »N 

A ( I  • J)  =  A ( I « J)  ♦  T*8(I*L) 

70  CONTINUE 

80  CONTINUE 

B<L*L>  =  -S*R 
DO  90  I =L 1 »N 
B(I.L)  =  0. 

90  CONTINUE 
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100  CONTINUE  ' 


reduce  a  to  upper  hessenberg.  keep  B  TRIANGULAR  ; 

IF  (N.LE.2)  GO  TO  170  l 

NM2=N-2  1 

DO  160  K*1,NH2  } 

K1  =  K*1  i 

NK1  =  N-K-l  3 

DO  150  LB=1*NK1  3 

L  =  N-L8  ] 

LI  =  L*1  \ 

CALL  HSR2(A(L.K) «A(L1.K).U1.U2«V1.V2)  3 

IF  (Ul.NF.l.)  GO  TO  125  j 

00  lio  J=K.N 

T  =  A (L. J)  ♦  U2*ACL1*J) 

A(L.J)  =  A (L« J)  ♦  T*V1  ] 

AILltJ)  =  A(L1* J)  ♦  T*V2  ] 

110  CONTINUE  ] 

A (LI .K)  s  0. 

DO  120  J=L.N  j 

T  =  B(L.J>  ♦  U2*B(LltJ) 

B (L* J)  =  B(L*J)  ♦  T*V1 

B (LI • J)  =  B(LltJ)  ♦  T*V2  j 

120  CONTINUE  ] 

125  CALL  HSH2(B(L1.L1).B(L1.L>,U1.U2.V1*V2> 

IF  (Ul.NE.l.)  GO  TO  150 
DO  130  1=1. LI 

T  =  B(I.Ll)  ♦  U2*B ( I *L)  ' 

B(I.Ll)  =  B(I.Ll)  ♦  T*V1 
B(I.L)  =  B(I.L)  ♦  T*V2 
130  CONTINUE 

6 (LI .L)  =  0. 

DO  140  1=1. N 

T  =  A(I. LI)  ♦  U2*A(I.L) 

A(I.L1)  =  Ad. LI)  ♦  T*V1 
A(I.L)  =  A(I.L)  ♦  T*V2 
140  CONTINUE 

IF  ( *NOT .WANTX)  GO  TO  150 
DO  145  1=1. N 

T  =  Xd.Ll)  ♦  U2*X(i.L) 

X(I.L1)  =  X(I.Ll)  ♦  T*V1 
X(i,L)  =  X ( I »L )  ♦  T*V2 
145  CONTINUE 

150  CONTINUE 
160  CONTINUE 
170  CONTINUE 

return 

ENO 
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SUBROUTINE  QZIT  (ND.N.A.B.EPS.EPSA.EPSB.ITER.WANTX.X) 
DIMENSION  A (NO »ND) *6 (NO. NO) .X(ND.ND) 

DIMENSION  ITER(N) 

LOGICAL  WANTX.MlD 
C 

C  INITIALI7F  ITER.  COMPUTE  EPSA.EPSB 
C 

ANORM  =  0. 
bNORM  =  0. 

DO  185  1=1. N 
ITEH(I)  =  0 
AN  I  =  0. 

IF  (I.NE.l)  ANI  =  ABS(AU.I-l) ) 

HNI  =  0. 

DO  180  J=I *N 

ANI  =  ANI  ♦  ABS(A( I . J) ) 

BN I  =  BN  I  ♦  ABS(B(I.J)) 

180  CONTINUE 

IF  (ANI. GT. ANORM)  ANORM  =  ANI 
IF  (BNl.GT. BNORM)  BNORM  =  BN I 
185  CONTINUE 

EPSA  =  EPS*AN0RM 
EPSB  =  EPS*BNORM 
C 

C  REDUCE  A  TO  QUASI -TRIANGULAR.  KEEP  B  TRIANGULAR 

C 

M  =  N 

200  IF  (M.LE.2)  GO  TO  390 


CHECK  FOR  CONVERGENCE  OR  REOUCIBILITY 

DO  220  LR=1.M 
L  =  M ♦ 1 “LB 

IF  (L.EQ.l)  GO  TO  260 

IF  (A6S(A(L«L-|))  .LE.  EPSA)  GO  TO  230 
220  CONTINUE 
230  A(L.L-l)  =  0. 

IF  (L.LT.M-1)  GO  TO  260 
M  =  L-l 
GO  TO  200 

CHECK  FOR  SMALL  TOP  OF  B 

260  IF  (ABS(B(L.L> > .GT.EPSB)  GO  TO  300 
B(L.L)  =  0. 

LI  =  L*1 

CALL  HSH2(A(L.L).A<L1.L).U1.U2*V1.V2) 

IF  (Ul.NE.l.)  GO  TO  280 
OC  270  J=L.N 

T  =  A(L.J)  ♦  U2*A (L  1 .  J) 

A(L.J)  =  A(L.J)  ♦  T*V1 
A(L1.J)  =  A(L1.J)  ♦  T*V2 
T  =  B (L. J)  ♦  U2*8(L1.J) 

B(L.J)  =  B(L.J)  «•  T*V1 
B(L1.J)  =  B(L1.J)  ♦  T*V2 
270  CONTINUE 
280  L  =  LI 

GO  TO  230 


36 


o  or-.  r  o  o  o  o  o  r, 


c 

C  BEGIN  ONE  QZ  STEP*  ITERATION  STRATEGY 

300  Ml  =  M  -  1 
LI  =  L  ♦  1 
CONST  =  0.75 
ITER (H)  =  ITER(M)  ♦  1 
IF  (ITER(M) .EO.l)  GO  TO  305 
IF  (A8S(A(M*M-1)).LT.C0NST*0LD1>  GO  TO  305 
IF  (ABS<A(H-1,M-2)).LT.C0NST*0LU2)  GO  TO  305 
IF  UTER(M).EQ.IO)  GO  TO  310 
IF  (ITER(M) .GT.30)  GO  TO  380 

ZEROTH  COLUMN  OF  A 

305  Bll  =  *»*•  .L) 

B22  =  B  <L1 tLl ) 

IF  (ABS(B22).LT.EPSB>  B22  =  EPSB 
B33  =  B(Ml.Ml) 

IF  (A8S(B33).LT.EPSB)  B33  =  EPSB 
844  =  B(M.M) 

IF  (ABS (844) .LT. EPSB)  B44  =  EPSB 
All  =  A (LfL)/Bl 1 
A12  =  A (L .Li )/B22 
A21  =  A(L1*L)/B11 
A22  =  A(L1*L1)/B22 
A33  =  A  <Mj  *Mj ) /B33 
A34  =  A (MJ  *M ) /B44 
A43  =  A (M.Ml )/B33 
A44  =  A (M.MJ/B44 
B 12  =  B(L,L1)/B22 
B34  =  8(M1*M)/B44 

A10  =  (  <A33-AU)*(A44-A11)  -  A34*A43  ♦  A43*B34*A11  )/A21 
I  ♦  A 1 2  -  A1 1*812 

A20  =  (A22-A11-A21*B12)  -  (A33-AH)  -  (A44-A11)  ♦  A43*B34 
A30  =  A (L*2*L 1 ) /B22 
GC  TO  315 

AD  hOC  SHIFT 

310  A10  =  0. 

A20  =  0. 

A30  =  1.1605 

315  OLOl  =  A8S(A(M*M-1 ) ) 

OLD2  =  ABS(A(M-l,M-2) ) 

IF  (.NOT.WANTX)  LORl  =  L 
IF  ( 4 ANTX )  LORI  =  1 
IF  (.NOT.rtANTX)  MORN  =  M 
IF  (WANTX)  MOkN  =  N 

BEGIN  MAIN  LOOP 

00  360  K-LtMl 
MIf)  =  K.NE.M1 
K1  =  K ♦ 1 
K2  =  Kt2 
K3  =  K*3 
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IF  (K3.GT.H)  K3  =  M 
KM1  =  K-l 

IF  (KM1.LT.L)  KM1  =  L 

IF  (K.EO.L)  CALL  HSH3(A10*A20,A30,U1 »U2»U3»V1 * V2, V3> 

IF  (K.GT .L, AND.K.LT.M1 ) 

CALL  HSH3 ( A (K »XM 1 ) , A  (K 1 *KM1 ) ,A (K2*KMl ) ,ul *U2»U3*  Vl » V2*  V3) 
IF  (K.£U.Hi>  CALL  HSH2 (A (KyKMj ) ,A CK j tKMl ) ,Uj »U2*Vl » Vp> 

IF  (Ul.NE.l.)  GO  TO  325 
DO  320  J=KM1 *MORN 

T  a  A(K»J)  ♦  U2«A(K1,J) 

IF  (MID)  T  =  T  ♦  U3*A(K2*J) 

A (K, J)  =  A(K.J)  ♦  T*V1 
A (K1 * J)  =  A ( K 1 1 J )  *  T»V2 
IF  (HID)  A(K2*J)  =  A(K2»J)  •»  T*V3 
T  =  B (K, J)  ♦  U2«B(K1«  J) 

IF  (MID)  T  =  T  ♦  U3*B(K2*J> 

B (K  * J)  =  B (K* J)  ♦  T»V  1 

B(K1, J)  =  B(K1,J)  ♦  T*V2 

IF  (MID)  B  (K2*  J)  =  B(K2»J)  *  T*V3 

continue 

IF  (K.EQ.L)  GO  TO  325 
A(K1,K-1)  a  0. 

IF  (MIO)  A(K2*K-1)  =  0. 

IF  (K.FU.M1)  GO  TO  340 

CALL  HSH3(B(K2»K2)*B(K2*Kl)'»B(K2*K)*Ul*U2*U3»VI»V2fV3) 

IF  (Ul.NE.l.)  GO  TO  340 
DO  330  I=L0R1,K3 

T  =  A (  I  *K2 )  ♦  U2*A ( J ,Kl )  ♦  U3*A(J,K> 

A(I,K2)  a  A(I,K2)  ♦  foVl 
A(I.Kl)  a  A(I.K1)  ,  T*V2 
A ( I *K)  =  A ( I *K )  ♦  T*V3 
I  =  0 ( I  *K2 i  ♦  U2*6 ( I  »Kl )  ♦  U3*B(I»K) 

B(I«K2)  =  B ( I *K2)  ♦  T*V1 
B(I.Kl)  a  8(1, Kl)  ♦  T«V2 
B(I.K)  a  B(I»K)  ♦  r*V3 
CONTINUE 

H(K2,K)  a  o. 

H(K2,K1)  a  0. 

IF  ( .NOT .WANT X)  GO  TO  340 
DO  335  1  =  1, M 

T  =  X ( I,K2)  ♦  U2#X ( I ,K1 )  ♦  U3*X ( I »K) 

X(I,K2)  =  X(I,K2)  ♦  T»V1 
X(I»K1)  =  X ( I ,K 1 )  ♦  T#V2 
X(I,K)  a  X ( J ,K )  ♦  T*V3 
CONTINUE 

CALL  KSH2 (B(K1 ,K1 ) , B ( K 1 ,  K  ) »U1 ,U2* Vl » V2 ) 

IF  (Ul.NE.l.)  GO  TO  36O 
DO  350  I=L0R1,K3 

T  =  A ( I ,K1 )  ♦  U2»A(I,K) 

A  ( I  ,K 1 )  =  A ( I . K 1 )  ♦  T«V1 
A  (  I  *K )  a  A  ( I  ,K )  ♦  T*V2 
T  =  8(1, Kl)  ,  U2*B ( I ,K ) 

B ( I ,K 1  )  a  B ( I  ,X1  )  ♦  T*v 1 
B(I,K)  =  8  ( I  , K )  ♦  T#V? 

CONTINUE 
R(K1,K)  a  0. 

IF  (.NOT.WANTX)  GO  TO  36(; 

DO  355  1=1, N 


non  non 


S®^5SSSS^VSS^2®y^ViWST.WS??WSJ,v-i SW7TV.> '.  v«^“!»’ 


T  =  X ( I «K  1 )  ♦  U2*X(I*K) 

X  ( I  *K  l )  =  X  ( I  »K  1 )  ♦  T*V1 
X ( I *K)  =  X(i» K)  ♦  T*V2 
355  CONTINUE 
360  CONTINUE 


.y; 


> 

„A 


END  MAIN  LOOP 


GO  TO  200 

END  ONE  QZ  STEP 

380  DC  385  1=1 «M 

ITEW(I)  =  -1 
385  CONTINUE 
390  CONTINUE 
pe  turn 
end 
c 
c 
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o  c.  o  non  r.  n  n  n  n  o  o  r> 


SUBROUTINE  QZVAL(ND*N,A*B*EPSA*EPSB*ALFR»ALFI»BETA*WANTX*X) 
DIMENSION  A(NDtND) *B(ND*NO) «ALFR(N) .ALFI (N) .BETA (N>  *X (ND«ND) 
LOGICAL  WANTX.FLIP 

FIND  EIGENVALUES  OF  QUASI -TRIANGULAR  MATRICES 

DO  400  THRU  490  FOR  M  =  N  STEP  <-l  OR  -2)  UNTIL  1 
M  =  N 

400  CONTINUE 

IF  (M.EO.l)  GO  TO  410 
IF  (A(M*M— 1) . NF • 0 . )  GO  TO  420 

0NE-8Y-0NE  SUBMATRIX*  ONE  REAl.  ROOT 

4 1 0  ALFK(M)  =  A (M*M) 

BETA(M)  =  B(M*M> 

ALFI(M)  =  0. 

M  =  M- 1 
GO  TO  490 

TWO-BY-TwO  SUBMATRIX 
420  L  =  M-l 

IF  (ABS(B(L*L) ) .GT.EPSB)  GO  TO  425 

8 (L*L)  =  0. 

CALL  HSH2 ( A (L *L ) *A(M*L>  *U1 *U2» VI  * V2> 

GO  TO  460 

425  IK  (ABS(B(M*M) ) .GT.EPSB)  GO  TO  430 

8  (M » M )  =  0. 

CALL  HSh2(A(M.M) »A(M*L) ♦U1*U2*V1*V2) 

8N  =  0. 

GO  TO  435 

430  AN  =  ABS(A(L*L))*A8S(A(L*M)) ♦ABS(A(M*L) ) *  ABS (A  <M*M) ) 

BN  =  ABS(B(L*L))*ABS<B<L»M>>*ABS(B(M*M>> 

Al)  =  A (L«L ) /AN 
A 1 2  =  A(L*M)/AN 
A21  =  A(M*L)/AN 
A22  =  A (M.M) /AN 
H 1 1  =  8(L*L)/BN 
d!2  =  B<L»V)/BN 
822  =  B  (M.M ) /BN 

C  =  (A11*B22  ♦  A22*B11  -  A21*B12>/2. 

D  =  ( A22*B 1 1  -  A1I*B22  -  A2 i *B1 2 > **2/4. 

1  ♦  A21*822*(A12*B11  -  A11*R12> 

IF  (D.LT.O.)  GO  TO  480 


TWO  REAL  ROOTS 

ZERO  BOTH  A (M*L)  AND  B(M,L> 

IF  (C.GE.O.)  E  =  (C  ♦  SQRT(D)  )/<Bll*B22) 

IF  (C.LT.O.)  E  =  (C  -  SORT (D) )/(Bll*B22> 

All  =  All  -  c«Rll 
A12  =  A 1 2  -  t*B12 
A22  -  A22  -  E*822 

FLIP  =  (ABS(All>*ARS(A12)).GE.(AbS(A21>*ABS(A22>) 
iF  (FLIP)  CALL  HSH2(A12»A11*U1*U2.V1*V2) 

IF  (.NOT. FLIP)  CALL  HSH2 ( A22*A21 *U1 ,U2*V1 *V2) 
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435  IF  (Ul.NE.l.)  GO  TO  450 
DO  440  1=1 ,M 

T  =  -.  ( 1  ♦  M )  »■  U2*A  ( I  »L ) 

A ( I , M )  =  A ( I .M )  ♦  V1*T 
A ( I «L )  =  A  (  I  tL )  ♦  V2*T 
T  =  B  <  I  *  M )  .  U2*B ( 1  «L ) 

B( I *M)  =  B ( I .M )  ♦  V 1 *T 
B { I ♦ L )  =  6 ( I «L )  ♦  V2*r 
440  CONTINUE 

IF  (.NOT.WANTX)  GO  TO  450 
DO  445  1=1. N 

T  =  X(I.M)  ♦  U2*X(1,L> 

X(I.M)  =  X ( I ?  M )  ♦  V 1*T 
X(I.L)  =  X(I.L)  ♦  V2*T 
445  CONTINUE 

450  IF  (HN.EO.O.)  GO  TO  475 

FLIP  =  AN  .GE.  A8S (E) *BN 

IF  (FLIP)  CALL  HSH2(B(L*L)*B(M*L)»Ul *U2* V 1  * V2 ) 

IF  (.NOT. FLIP)  CALL  HSH2 ( A (L.L) ♦ A (H.L) *U1 .U2.V1 .V2) 

460  IF  (Ul.NE.l.)  GO  TO  475 

DO  470  J=L«N 

T  =  A(L.J)  ♦  U2»A(M.J) 

A(L.J)  =  A (L ♦ J)  ♦  V1*T 
A(M,j)  =  A(M.J)  ♦  V2*T 
T  =  B(L.J)  ♦  U2«8(M,J) 

B(L.J)  =  B(L.J)  ♦  V 1 *T 

8(M,J)  =  B(M.J)  ♦  V2*T 

470  CONTINUE 

475  A ( M . L )  =  0. 

h(M.L)  =  0. 

ALFB(L)  =  A(L.L) 

ALFP(H)  =  A (M.M) 

BETA(L)  =  B (L.L) 

BETA(M)  =  B(M,M> 

ALFI(M)  =  0. 

ALFI  (L)  =  0. 

M  =  M-p 
GO  TO  490 
C 

C  Two  COMPLEX  HOOTS 

C 

4« 0  E*  =  C/  ( 61  1  ^2? ) 

El  =  SGRT(-D)/(Bll«6?2> 

A I  1 H  =  AH  -  ER*Bll 
a l l I  =  EI*P11 
A12H  =  A 1 2  -  ER*H12 
A12I  =  EI*612 
A  2 1 h  =  A21 
A? 1  I  =  0. 

A22R  =  A22  -  E H * B 2 2 
“221  =  Ei*H22 

FLIP  =  (ABS(A11R)*ABS(A1U>*ABS(A12RMABS(A12I>>  .Gt. 

1  (ABS (A21R) ♦A8S(A2?K)»ABS(A22I) ) 

IF  (FLIP)  CALL  CHSP2(A12R.A12I.-A11R.-A11I.C2.SZR.S2I) 

IF  (  .NOT  <,FL  I P )  CALL  CHSH? ( A22R ♦ A22 1 ,-A2 1R *-A21 1 .CZ ♦ SZR t SZ I  - 
FLIP  =  An  .GE.  (ABS (EH) ♦ABS(El) )»BN 
IF  (FLIP)  CALL  CHSH2(CZ»H11*SZR*812,  SZl*B12. 

1  SZR*B22»  SZI*R22.  CG,  SOR.  SQI) 
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JPSEfS 


IF  (.NOT. FLIP)  CALL  CHSH? (CZ*A1 1 ♦SZR*A12»  SZI*A12* 
1  CZ*A21 ♦SZR*A22*  SZI*A22*  CQ*  SQR*  SQI) 

SSR  =  SQR*SZ«  ♦  SQI*SZI 
SSI  =  SQR*SZI  -  SQI*SZR 

TR  =  CQ*CZ*Al 1  ♦  CQ*SZR#Ai2  ♦  SQR*CZ#A21  ♦  SSR*A22 
T I  =  CQ*SZI*A12  -  SQI*CZ*A21  ♦  SSI*A22 
8DR  =  CQ*CZ*B11  ♦  CQ*SZR*812  ♦  SSR*B22 
BDI  =  CQ*SZI«B12  ♦  SSI*B22 
R  =  SORT (BDR*BOR  ♦  BDI*BDI) 

BETA (L )  =  BN«R 

ALFR(L)  =  AN»(TR*B0R  ♦  TI*BDI)/R 


ALT  a:'.)  =  AN*(TR*8DI  -  TI*BDR)/R 

TR  =  SSR*A1I  -  SQR*CZ*A1?  -  CQ*SZR*A21  ♦  CQ*CZ*A22 
T I  =  SSI*A11  -  SQI*CZ*A12  ♦  C0*SZI*A21 
BOR  =  SSR*B11  -  SQR*CZ*B12  ♦  CQ*CZ*B22 
BDI  =  -  SS1*B11  -  SQI*CZ*B12 

R  =  SORT (BDR*80R  ♦  BDI*BDI ) 

BETA(M)  =  BN*R 

ALFR(M)  =  ANMTR*BDR  ♦  TI*BDI)/R 
ALFI(M)  =  A*-'"  <TR*BDI  -  TI*BOR)/R 
M  =  H -2 

C 

<♦90  IF  (M.GT.O)  GC  TO  400 
RETURN 
E.%0 


C 


ooo  n n o n  n n 


ALFM  =  ALFR(M) 

BETM  =  BETA (M) 

IF  (ABS(ALFH) .LT.EPSA)  ALFM  =  0. 

IF  (aBS(BETM) .LTvEPSB)  BETH  =  0. 

B (M.M)  -  1. 

C 

C  DO  510  THRO  540  FOR  L  =  M-l  STEP  (-1  OR  -2>  UNTIL  1 

C 

L  =  M-l 

IF  (L.EQ.O)  GO  TO  540 
510  CONTINUE 

LI  =  L ♦ 1 
SL  =  0. 

DO  515  J=L1 «M 

SL  =  SL  ♦  (BETM*A(L»J>-ALFM*B(L*J>)*B(J»M) 

515  CONTINUE 

IF  (L.EQ.l)  GO  TO  520 
IF  (A(L»L-1) .NE.O.)  GO  TO  530 
520  D  =  BETM*A(L*L) -ALFM*B (Lt  L) 

IF  (D.EQ.O.)  D  =  (EPSA*EPSB)/2. 

B(L*M)  =  -SL/O 
L  =  L-] 

GO  TO  540 
C 

530  K  =  L- 1 

SK  =  0. 

00  535  J=L i « M 

SK  =  SK  ♦  <BETM«A(K»J)-ALFM*B(K*J) >*B(J»M> 

535  CONTINUE 

TKK  =  BE  TM*A  (K  »K  )  -  Al_FM*8(K,K) 

Tru  =  eETM»A(K,L)  -  ALFM«B(KfL) 

TLK  =  BETM»A(L»K) 

TLL  =  BETM*A(L»L)  -  ALFM*B (L.L) 

0  =  TKK*TLL  -  TKL^TLK 
IF  (O.EC.C.)  0  =  (EPSA*EPSB) /?• 

B(LsM)  =  (TLK*SK  -  TKK*SL)/D 

flip  =  absitkki  .ge.  absitlkj 

IF  (FLIP)  8(K,M)  =  -(SK  ♦  TKL*B(LtM) )/TKK 
IF  (.NOT. FLIP)  B(K.M)  =  -(SL  ♦  TLL*B (L »M> ) /TLK 
L  =  L-? 

540  IF  (L.GT.O)  GO  TO  510 
M  =  M-l 
GO  TO  590 
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COMPLEX  \/ EC  I  OR 

ALMR  =  ALFR(M-l) 

ALMI  =  ALFI (M-l) 

BE  TM  =  BtTA(M-l) 

MR  =  M-l 
MI  =  M 

NORMALIZE  SO  THAT  M-TH  COMPONENT  =  0.~1«*I 

(M-l)ST  =  -<BETM*A<M,M>-ALFM*B<M,M)>*(M-YH>/<BETM*A(M,M-1>  > 

H (M-l «  MR )  =  ALMI*B(M,M)/ibEfM*A(M*M-l)) 

B(M-1,M1)  =  (BETM*A(M«M)-ALMR«B(M*M) )/(BETM*A(M*M-l ) ) 

M(M»MR)  =  0. 

M (MfMI )  =  -1. 

DO  560  Thru  585  FOR  L  =  M-2  STEP  (“1  OR  -2)  UNTIL  1 


IF  (L.EQ.O)  60  TO  585 
CONTINUE 
LI  =  L*1 
SLR  =  0. 

SLI  =  0. 

DO  565  J=L1,M 

TR  =  BETM*A(L.J)  -  ALMR*B(L,J> 

T:  =  -ALMI«B(L*J) 

SLR  =  SLR  ♦  TR*B( J*MR)  -  TI*B(JtMI> 

SLI  =  SLI  ♦  TR*B ( J»MI )  ♦  TI*B(J.MR) 

CONTINUE 

IF  (L.EO.l)  GO  TO  570 
IF  (A(L.L-l) •NE.O.)  GO  TO  575 
DR  =  BETM»A(L*L)  -  ALMR*B(L*L) 

01  =  -ALM 1*8 (L  tL ) 

CALL  CDIV(-SLR,  -SLI,  DR»  Dl*  B(L*MR)  ♦  B(L,Ml)> 
L  =  L-l 
00  TO  585 

K  =  L-l 
SKk  =  0. 

SKI  =  0. 

DO  580  J=L l iM 

TR  =  B£TM*A(KtJ)  -  ALMR*B<K,J> 

TI  =  -ALMI*B(K,J) 

SKR  =  SKR  ♦  TR*8 ( J»MR )  -  TI*B(JfMl) 

SKI  =  SKI  ♦  TP*B(J,KI>  ♦  T I *B ( JtMR) 

CONTINUE 

TKKR  =  BETM*A(K,K)  -  ALMR*B (K *K ) 

TKKI  =  -ALMI*B(K»K) 

TKLR  =  BETM*A (K  «L )  -  ALMR«B(K,L> 

TKL I  =  -ALM I *B  <K  »L ) 

TLKR  =  BETM*A(L,K) 

TLK I  =  0. 

TLLR  =  8ETM*A(L*L)  -  ALMR*B(LtL) 

TLLI  =  -ALM I *B (L ,L ) 

DR  =  TKKR*TLLR  -  TKKI*TLLI  -  TKLR*TLKk 
01  =  TKKR*TLLI  ♦  TKKI*TLLH  -  TKLI*TLKR 


n  r.r  nor-. 


^vmsnmm 


Vmy: 


1 

2 


1 

2 

1 

2 


585 

590  IF 


IF(DR.EG.0.  .AND.  0I.EQ.0.)  DR  =  <EPSA*EPSB)/2. 

CALL  CDIV(TLKR*SKR-TKKR*SLR*TKKI*SLI. 

tlkr*ski-tkkr«sli-tkki*slr, 

OR,  DI,  B(L,MR),  B (L,MI ) ) 

FLIP  =  (ABS (TKKR) ♦  ABS (TKKI ) )  .GE.  ABS(TLKR) 

IF  (FLIP)  CALL  CDIV(-SKR-TKLR*B(L,MR)>TKLI*8(L*M1) , 

-SKI-TKLR*B(L»MI)-TKLI*B(L,MR) . 

TKKR,  TKKI,  B(K,MR),  B(K»MI)) 

IF  (.NOT.  FLIP)  CALL  CD IV (-SLR-TLLR*B (L »MR) ♦TLLI*B(L»MI) 

-SLI-TLLR*B(L,MI)-TLLI*B(L,Mk)  , 

TLKR,  TLKI,  B(K,MR),  8  (K,MI ) ) 

L  =  L-p 

IF  (L.GT.O)  GO  TO  560 
M  =  M-2 

(M.GT.O)  60  TO  500 


TRANSFORM  TO  ORIGINAL  COORDINATE  STSTEM 
M  =  N 

600  CONTINUE 

00  6?0  1=1, N 
S  =  0. 

00  610  J=1,M 

S  =  S  ♦  X ( 1 « J)*B( J«M) 

610  CONTINUE 

X  ( I  ,  M )  =  S 

620  CONTINUE 
M  =  M-l 

IF  (M.GT.O)  GO  TO  600 
NORMALIZE  SO  THAT  LARGEST  COMPONENT  =  1. 


M  =  N 

630  CONTINUE 
S  =  0. 

IF  (ALFI (M) .NE.O.)  GO  TO  650 
00  635  1=1. N 

R  =  ABS ( X ( I »M ) ) 

IF  (R.LT.S)  GO  TO  635 
S  =  R 

0  =  X  ( I ,  M ) 

635  CONTINUE 

00  640  1  =  1 »  N 

X ( I ,M)  =  X(I,M)/0 
640  continue 

M  =  M-l 

GO  TO  690 

650  00  655  1  =  1. N 

R  =  X ( I ,M-1 ) **2  ♦  X ( I ,M ) n#2 
IF  (R.LT.S)  GO  TO  6cm 
S  =  R 

OR  =  X(I, M-l) 

01  =  X  ( I , M ) 

655  CONTINUE 

00  660  1=1, N 

CALL  CDIV(X(I,M-l),X(l,M),DR,DI,Xil,M-l),XfI*M>) 
660  CONTINUE 


45 


iwpwiSiWiisws’iJ  .r- 


.-.  v  ^VT.  AV.  ^-'/ -'/v-^1' *! ''!  -.--  ■-  -  ■  v'-T-  ,^-’C2\335reg,JE5Sra 


M  =  M-2 

690  IF  (M.6T.O)  GO  TO  630 

C 

700  RETURN 
ENO 
C 

c 


] 

i 
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SUBROUTINE  HSH3<Ai,A?,A3*U1»U2»U3*V1*V?*V3) 

FINOS  HOUSEHOLDER  TRANSFORMATION  THAT  WILL  ZERO  A?  AND  A3 
P  =  I  ♦  (Vl,V2»V3)*(Ul*U2tU3)*«T 

IF  (A2  EG.Oc  .AND.  A3.EQ.0.)  GO  TO  10 
S  =  ABS(Al)  ♦  ABSIA2)  ♦  ABS(A3) 

U1  =  Al/S 
U2  =  A2/S 
U3  =  A3/S 

R  =  SORT (U1*U1 ♦U2*U2^U3*U3> 

IF  (U1.LT.0.)  R  =  -R 
VI  =  -<U1  ♦  R)A 
V2  =  -U2/R 
V  3  =  -U3/R 
U1  =  1. 

U2  =  V2/V1 
U3  =  V3/V1 
RtTuRN 
10  U1  =  0. 

RETURN 

END 


i 
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SUBROUTINE  HSH2< Al *A2«Ul *U2»V1 *V2) 

FINOS  HOUSEHOLDER  TRANSFORMATION  THAT  WILL  ZERO  A2 
p  =  I  ♦  <V1,V2)*(U1*U2)**T 

IF  (A2.EQ.0.)  GO  TO  10 
S  =  ABS(Al)  ♦  A8S(A2) 

U1  =  Al/S 
U2  =  A2/S 

R  =  SORT <U1*U1*U2«U2) 

IF  (U1.LT.0.)  R  =  -R 
VI  =  -<U1  ♦  R)/R 
V2  =  -U2/R 
U1  =  l. 

U2  =  V2/V1 

Return 
lo  ui  =  o. 
return 
End 
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SUBROUTINE  CHS M2 < AIR* Ai I  * A2R* A2I *C*SR*  SI ) 

COMPLEX  HOUSEHOLDER  THAT  WILL  2EHO  A? 

(C  S*> 

P  =  (S  -C  )  *  C  REAL*  S  COMPLEX 

IP  (A2R.EQ.0.  « AND •  A2I.EQ.0.)  GO  JO  10 
IP  (AlR.EQ.O.  .AND.  A1I.EQ.0.)  GO  TO  20 
R  =  SQRT(A1R*A1R^A1I*A1I) 

C  =  R 

SR  =  (AlR*A2R*AlI*A2D/R 
SI  =  (AiR*A2l-AlI*A2R)/R 
R  =  SORT<C*C*SR»SR*SI*Sl) 

C  =  C/R 
SR  =  SR/R 
SI  =  SI/R 
Rfc TURN 
10  C  =  1. 

SR  =  0. 

SI  =  0. 

RETURN 
20  C  =  0. 

SR  =  l. 

SI  =  o. 

RETURN 
END 

C 

C 


49 


AoaaiK: 


non 


SUBROUTINE  CDIV(XRtXl«YR»Yl*ZR*ZI > 


COMPLEX  DIVIDE*  Z  =  X/Y 

IF  (ABS(YR) .LT.ABS(YI) >  GO  TO  10 
MR  =  XR/yR 
WI  =  XI/YR 
VI  =  Yl/YR 

D  =  1.  ♦  VI*VI 
ZH  =  <WR  ♦  WI*V1)/D 
Z I  =  (W I  -  WR*VI)/D 
RETURN 

10  MR  =  XR/YI 
Ml  =  XI/YI 
VR  =  YR/YI 
D  =  VR*VR  ♦  1. 

ZR  =  (WR*VR  ♦  WI)/D 
ZI  =  <WI*VR  -  WR)/D 
return 

END 


